Simulation of fluid-solid coexistence in flnite volumes: A method to study the 
properties of wall-attached crystalline nuclei 
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The Asakura-Oosawa model for colloid-polymer mixtures is studied by Monte Carlo simulations 
at densities inside the two-phase coexistence region of fluid and solid. Choosing a geometry where 
the system is confined between two flat walls, and a wall-colloid potential that leads to incomplete 
wetting of the crystal at the wall, conditions can be created where a single nanoscopic wall-attached 
crystalline cluster coexists with fluid in the remainder of the simulation box. Following related 
ideas that have been useful to study heterogeneous nucleation of liquid droplets at the vapor-liquid 
coexistence, we estimate the contact angles from observations of the crystalline clusters in thermal 
equilibrium. We find fair agreement with a prediction based on Young's equation, using estimates 
of interface and wall tension from the study of flat surfaces. It is shown that the pressure versus 
density curve of the finite system exhibits a loop, but the pressure maximum signifies the "droplet 
evaporation-condensation" transition and thus has nothing in common with a van der Waals-like 
loop. Preparing systems where the packing fraction is deep inside the two-phase coexistence region, 
the system spontaneously forms a "slab state" , with two wall-attached crystalline domains separated 
by (flat) interfaces from liquid in full equilibrium with the crystal in between; analysis of such states 
allows a precise estimation of the bulk equilibrium properties at phase coexistence. 
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I. INTRODUCTION AND OVERVIEW 

Nucleation of crystals from fluid phases is an impor- 
tant problem [1-5 with important applications, such as 
formation of ice crystals in the atmosphere, solidification 
of molten silicates in processes deep underneath the earth 
crust, and last but not least crystallization processes of 
various materials are ubiquitous in many technical pro- 
cesses. However, nevertheless crystal nucleation is rather 
poorly understood on a quantitative level: one mostly re- 
lies on the concept of classical nucleation theory [2|[6UTQ]. 
but since almost always the "critical nucleus" that trig- 
gers the phase transition contains only a few tens to at 
most a few thousand particles, considerations based on 
macroscopic concepts (balancing bulk and surface free 
energies, using the interfacial tension of macroscopic flat 
interfaces, etc. [1-5 ) are doubtful. Moreover, in most 
cases of interest nucleation is not homogeneous (i.e., trig- 
gered by spontaneous thermal fluctuations) but rather 
heterogeneous [9l [TIHT5] (i.e., triggered by defects, e.g. 
nucleation of "droplets" attached to a flat solid wall fa- 
cilitating formation of crystal planes parallel to the wall 
stacking on top of each other). Again, macroscopic con- 
cepts (involving the "contact angle" pIBUTS] at the wall, 
and possibly a free energy excess due to the three-phase 
contact line where the crystal and the fluid meet at the 
wall, the "line tension" [16, 19 - 21 are used [22 , but their 
reliability is uncertain. 

Recently progress has been achieved by studying the 
nucleation of colloidal crystals [23-27 . Colloidal par- 
ticles are in the jim range, and hence the structure of 
fluid-crystal interfaces can be studied with single-particle 
resolution [28 , and since dynamics of such systems are 
very slow, the time evolution of interfacial phenomena 
can be followed in real time |29l |30| • A further bonus is 



that effective interactions between colloidal particles are 
tunable to a large extent J3TH33] . A good example of this 
point are colloid-polymer mixtures [34-36 : varying the 
polymer concentration in the colloidal dispersion one can 
vary the depletion attraction between the particles [^ . 

This effect has flrst been described in term of a simple 
model, the Asakura-Oosawa model J37H39] . and subse- 
quently it has been shown that this model does account 
qualitatively for the experimental observations very well 
|35l [36^. Thus, varying the polymer concentration the 
width of the two-phase coexistence region and the mag- 
nitude of the interfacial tension between coexisting fluid 
and solid phases can be controlled [40j, as sketched in 
Fig. [T] If the colloid-polymer mixture is conflned be- 
tween two (equivalent) walls, and the walls are prepared 
such that there is incomplete wetting [16-18 of the solid 
at the walls, it is likely that by variation of the polymer 
concentration in the dispersion one can change the con- 
tact angle. Such a control of wetting properties is very 
difficult to achieve in small molecule systems. 

In the present paper we hence want to contribute to 
the theoretical understanding of wall-attached crystalline 
clusters in colloid-polymer mixtures by computer simula- 
tion methods. As is well known, for macroscopic systems 
the "critical droplets" occurring in nucleation processes 
are hard to observe, since one has to focus on transient 
rare events when the system traverses a saddle point in 
the free energy landscape [4TH43] . when a droplet grows 
from subcritical to supercritical size. Here, computer 
simulations possess an advantage because in a system 
of conserved density in a finite-sized simulation box the 
coexistence of the critical droplet with surrounding "par- 
ent" phase is a situation of stable equilibrium [44J, un- 
like the situation in the thermodynamic limit where the 
"parent" phase is metastable, and the droplet on top of 
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FIG. 1. Schematic phase diagram of the Asakura-Oosawa 
model for a colloid polymer mixture in the plane of variables 
polymer reservoir packing fraction rf!p and colloid packing frac- 
tion Tj (left part), and variation of the fluid-solid interface 
tension jioocr'^ /ksT vs. rjp (right part). Here a is the colloid 
diameter, ap the polymer diameter, q = dp /a and q^ — 0.154 
the critical value of g, above which three-body interactions 
appear in the effective interaction between colloids mediated 
by the polymers. The region of two-phase coexistence, in be- 
tween the volume fraction n/ (where freezing begins) and rfrn 
(where melting begins) has been shaded. Note that for fluid- 
solid interfaces the interface tension depends on the interface 
orientation, 7100 refers to an interface that is perpendicular 
to the X-axis. 



the saddle is in unstable equilibrium 03]. In previous 
work, it has been shown that critical droplets (or bub- 
bles, respectively) associated with the liquid-vapor tran- 
sition [45^ ^ or systems undergoing an unmixing tran- 
sition in symmetric binary fluids [46l [47] or Ising models 
[HI [HI [48] can be studied in this way. When the pack- 
ing fraction {rf) of the finite system is chosen such that 
it falls inside of the two-phase coexistence region of the 
infinite system, i.e. rjf < r] < r]m (Fig. IT]), we may en- 
counter phase coexistence inside the simulation box. For 
a fluid-solid transition in thin film geometry with walls 
where incomplete wetting by the crystal occurs, we ex- 
pect various shapes of the minority domain, depending 
on T] as shown in Fig. |2J Of course. Fig. |2] is inspired 
by analogous studies of vapor-liquid transitions [HI [T^ , 
where both coexisting fluid phases are homogeneous and 
isotropic: only then it does clearly make sense to de- 
scribe the minority domain as sphere caps (Fig. [2k) or 
cylinder caps (Fig. ^p) that are stabilized by the peri- 
odic boundary condition (and oriented in the x-direction 
when Lx < Ly and along the ^/-direction when L^ > Ly^ 
while for L^ = Ly there occurs a degeneracy). For solid- 
liquid coexistence in finite volumes, Fig.[2]is approximate 
because of two reasons: (i) on the nanoscale, when the 
height of the sphere cap or cylinder cap is only a few lat- 
tice spacings, the discrete lattice structure of the crystal 
should be taken into consideration (ii) only when the lin- 
ear dimensions of the crystalline domain are very much 
larger than the lattice spacing and therefore the question 
what is its "macroscopic shape" makes sense. But even 
then Figs.[2|i), b) only hold when the fluid-solid interface 
tension does not depend on the orientation of the inter- 
face. Already in the bulk the shape of a crystal in general 
therefore is never a sphere, but rather needs to be found 



from the anisotropic interface tension via the Wulfl con- 
struction ^49. - .51j . The extension of this construction to 
wall-attached crystals has been given by Winterbottom 
et al. [52'-'5T. In fact, it is a nontrivial question under 
which conditions planar facets (rather than curved in- 
terfaces) occur [55 . Of course, for nano-crystallites we 
also expect a finite-size rounding of faceting transitions, 
related to the finite-size rounding of the interfacial rough- 
ening transition [56] , and hence the analysis of the equi- 
librium shapes of nano-crystals is very subtle. In Fig. [2| 
we have also assumed that the geometry (and conditions 
at the wall that occurs dit z = D) can be chosen such 
that no interfaces occur that connect both walls (though 
in fluid systems the occurrence of "liquid bridges" is ex- 
tremely common, see e.g. [57-61 ). Thus, the present 
study can present first exploratory steps only. 
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FIG. 2. Schematic description of two-phase configurations 
in a thin film geometry of volume Lx x Ly x D geometry, with 
periodic boundary conditions in x- and ^-directions, and two 
walls at z = and z — D. It is assumed that the wall at z = 
exhibits incomplete wetting by the crystal, while the wall at 
z — D IS not shown. All three cases shown assume a colloid 
packing fraction 77 in the two phase coexistence region chosen 
such that the solid phase is the minority phase. For rjt < rj < 
rjt' the solid domain is a sessile sphere-cap-shaped crystalline 
cluster attached to the wall (a) while for r]t' < r/ < r/t" it 
is a cylinder cap (b), and for 77 > 77^// it is a crystalline slab 
(c). Note that cases (a,b) ignore the crystalline structure of 
the solid, which rather is treated in a continuum description, 
ignoring both the anisotropy and discreteness of the lattice 
structure 

In Sec. II, we shall define precisely the model that 
is simulated and give details on the techniques of "sys- 
tem preparation" , simulation methods, and analysis tech- 
niques (note that it is a delicate matter to decide which 
particles are to be counted as part of the crystalline solid 
domain or as part of the fluid in each microstate of the 
simulation [40 ). In Sec. Ill we describe our numerical 
results and discuss them in the light of the questions that 
have been outlined above, while Sec. IV summarizes our 
conclusions. In the Appendix A, we consider the coexis- 
tence between wall-attached crystalline films and a fluid 
state in between, separated from the crystalline layers 
by planar interfaces. This is a "self-regulating" system. 



where the thickness of the fluid phase adjusts itself, so 
that the lever rule holds. We show that this geometry 
is useful for a direct estimation of the bulk packing frac- 
tions 77/, r]m and the coexistence pressure Pco- In the ap- 
pendix B, we compare the metastable crystallites (with 
100 faces adjacent to the wall) to stable ones (where the 
close-packed 111 faces are adjacent to the wall). 



II. MODEL, SIMULATION AND ANALYSIS 
TECHNIQUES 

The Asakura-Oosawa model of colloid-polymer mix- 
tures [37H39] describes the colloids as hard spheres of 
diameter a, the polymers are described as soft spheres of 
diameter ap^ and both colloid-colloid and colloid-polymer 
overlap is strictly forbidden, while polymer-polymer over- 
lap does not cost any energy. Thus, the interaction po- 
tentials are ("c" stands for colloids, "p" for polymers, r 
is the distance between the particles) 



Ucc{r < cr) = 00 , Ucc{r > a) =0 



(1) 



Ucp{r < (a + ap)/2) = 00 , Ucp{r > (a + ap)/2) = , 



and 



Upp{r) = 



(2) 



(3) 



The packing fraction of polymers {r]p) and colloids {r]) 
are then deflned in terms of the corresponding densities 
pp = Np/V and pc = Nc/V of these particles (V is the to- 
tal volume, A/p, Nc are the particle numbers of polymers 
and colloids, respectively) 



r]p = (7rcr^/6)pp , T] = {7ra^/6)pc 



(4) 



It is convenient to use the chemical potential /ip of the 
polymers as an external control variable, or, equivalently, 
the fugacity Zp = exp(/ip//c^T). The "polymer reservoir 
packing fraction" rj^ then is deflned as 



Vp = {^^p/^)zp 



(5) 



Being interested in static equilibrium properties of the 
model, one can proceed by flrst integrating out all the 
coordinates of the polymers, keeping only the coordinates 
of the colloids in the system as variables. For q = Opja < 
q^ = 0.154, this can be done explicitly and shown to yield 
an effective colloid-colloid attraction fMl ESl 



a < r < a + ar, 



3r 



2o-(l + q) 



2cr3(l + g)3 



(6) 



Ucc{r > cr + (Tp) = 



(7) 



Of course, for r < a we still have Eq. (IT]); Eqs. (p|, (ItI 
fully account for the depletion attraction between the col- 
loids caused by the polymers and show that the strength 
of this interaction can easily be controlled by variation 
of rjp. In the present study, we choose a single value of 
r]p and a single choice of q only. 



Vv 



0.1, q = 0.15 



(8) 



We also note that for 77^ = the model reduces to 
the simple hard sphere model, for which interfacial prop- 
erties (see e.g. [40l [64] for references) and nucleation 
(e.g. [T2I [13^) have been studied extensively; but since 
there is evidence [65] that hard spheres at hard walls 
(as well as on walls where a soft repulsion acts [64 ) ex- 
hibit complete wetting when the freezing fraction rjf (cf. 
Fig. [1]) is approached, the simple hard sphere model is 
less suitable to study crystalline nuclei attached to flat 
walls, and shall not be considered here further. 

As already indicated in Fig. [2] we choose a LxXLyX D 
geometry with periodic boundary conditions in x and 
y directions, while soft repulsive walls occur at z = 
and z = D^ respectively. These walls are described by a 
potential of the Weeks- Chandler- Andersen \^6^ type 



{kBT)-'VwcA{z) = 4e[{a^/zy^ 

for < ^ < cr^2^/^ , 

= 4e[{a^/{D - z))'' - {a^/{D - z)f 

for(D-cr^2^/^) <z<D , 

= otherwise. 



{a^/zf + 1] 

1] 



(9) 



Here e describes the strength of the potential (in units 
of the thermal energy ksT) and cr^ its range; in the 
present paper we only consider the case e = 1^ a^ = a 12. 
Choosing a = 1 as our unit of length, typical box linear 
dimensions were 



Lx = Ly= 39.606548352, D = 40.313808144 (10) 

which means that for a typical packing fraction rj = 
0.511184 those linear dimensions correspond to a total 
number of colloids Nc = 61717 in the system (allowing 
the observation of crystalline "clusters" containing 4743 
colloids, for instance, see below). 

In addition to simulation boxes where the wall surface 
has square geometry, we have also chosen L^ , Ly consis- 
tent with a perfect triangular lattice (of the lattice spac- 
ing corresponding to r]m)^ to allow the formation of fee 
crystalline layers with close-packed planes at the wall. 

An important ingredient in our study is the accurate 
knowledge about the phase transition in the bulk for our 




t being the time variable in the simulation: if p < Pco, 
the crystal shrinks on average (and V increases), while 
for p > Pco the crystal grows. This method has been care- 
fully tested for the simple hard sphere model [40", "W, and 
found to give very accurate results; we expect it to work 
for the AO model similarly well. 
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FIG. 3. Pressure p (in units of ksT/a^) plotted vs. pack- 
ing fraction, from the AO model with q = 0.15, r]p = 0.1. 
This phase diagram was already obtained by Zykova-Timan 
et al. [40], who found that fluid-solid coexistence occurs at 
Pco = 8.00 (highlighted by the dotted line). For 77 slightly 
larger than 77/, data points below the metastable bulk fluid 
branch (shown in the insert on largely magnified scales) corre- 
spond to states of two-phase coexistence, where a solid wall- 
attached crystalline cluster coexists with fluid phase, as shown 
schematically in Fig. [2^. Note that 77/ = 0.494 and ?7m = 0.64, 
respectively. 



model. We have made the particular choice of Eq. ([8| 
because for this choice the phase transition has already 
been studied by Zykova-Timan et al. [40j using constant 
pressure Monte Carlo methods (NpT ensemble [67| l68]). 
Fig.[3]shows the resulting equation of state: there are two 
branches of the pressure versus packing fraction curve, 
a fluid branch and a solid branch, which corresponds to 
the face-centered cubic (fee) lattice structure. These data 
were obtained from simulations of cubic LxXLyX L^ sim- 
ulation boxes with periodic boundary conditions, with 
Lx = Ly = Lz = L for the fluid branch, while the ratios 
of Ly/Lx and L^/Lx for the simulation of the crystal were 
adjusted such that an integer number of close-packed 
lattice planes stacked upon each other in the fee AB- 
CABC stacking sequence were compatible with the peri- 
odic boundary conditions without distorting the lattice. 
As is evident from Fig. [3J there occurs hysteresis over a 
broad range of pressures. The estimation of the pressure 
Pco at which phase coexistence occurs in equilibrium was 
done [40] using a method described in Refs. [40l[69]. In 
short, in this method one prepares a slab configuration, 
where in an elongated simulation box {L x L x 5L) a 
crystalline domain is separated from liquid domains to 
the right and to the left by flat domain walls running 
perpendicular to the z-direction. Precautions are taken 
to avoid any elastic deformation of the crystal slab in 
such simulations. Then the average volume of the system 
(V(t)) is determined in an NP^T Monte Carlo simulation 
as a function of Monte Carlo time for various pressures, 
and Pco is found from the condition that d{V{t))/dt = 0, 



FIG. 4. Initialization of the system: To realize a single crys- 
talline cluster attached to the wall at z = 0, we prepare the 
total system first at packing fraction 77/, and then replace ei- 
ther a cuboid (left) or a hemisphere (right) by a face centered 
cubic crystal, with the close packed (111) planes stacked par- 
allel to the plane z = 0. The volume of the cuboid (or hemi- 
sphere) is chosen in accord with Eq. ( 13 ) at each chosen value 
of 77. 



As a complementary approach to obtain the equation 
of state from simulations in the NpT ensemble, we have 
developed a method to very accurately compute the pres- 
sure from simulations in the constant volume (NVT) en- 
semble [70 . While the contribution to the pressure from 
the attractive part of the potential {Eq. (|6|} is straight- 
forwardly obtained from the standard virial expression 
j67l [68] , care is needed for the accurate estimation of 
the pressure contribution due to the hard core repul- 
sion, Eq. ([l) [70 . In a recent paper, we have shown 
that this task can be solved by applying the method due 
to DeMiguel and Jackson [TlJ. In this way one can ob- 
tain both the hydrostatic pressure of a bulk system as a 
function of the packing fraction, and for a system with 
external walls one can extract the wall tension as usual 
[T6] from the anisotropy of the pressure tensor. We re- 
fer the interested reader to Refs. [6T, TO for details on 
this method. Here, we only emphasize the two following 
facts: (i) for a homogeneous bulk system, the function 
p{r]) found in the NVT ensemble precisely coincides with 
its counterpart in the NpT ensemble [70 : thus finite size 
effects due to the change of the statistical ensemble are 
negligibly small in our problem, (ii) The method of Ref. 
[70] yields also the local transverse component Pt{z) if 
we deal with a system confined by walls. This allows to 
estimate the pressure Pf{r]) in the fluid part of a system 
that has separated into a solid cluster and surrounding 
fluid. Since this fluid coexisting with a crystalline clus- 
ter of finite size does not have the coexistence pressure 
Pco of the bulk, but rather the pressure of the fluid co- 
existing with the crystalline cluster is enhanced (Laplace 
pressure), estimation of this pressure is nontrivial and of 
interest: these data hence are included in Fig. [3J but we 
defer their analysis to the next section. 



At this point, we emphasize that it is rather straight- 
forward to prepare (Appendix A) a system at packing 
fraction near (77^ + 77^)/2, the center of the two-phase 
coexistence region, then the system develops easily to- 
wards an equilibrium state, where both walls are coated 
by thin crystalline films (Fig. ^) with a liquid state of 
appropriate thickness in between, and the pressure in- 
side the liquid must be the coexistence pressure Pco- By 
this method, the results of [40 could be checked inde- 
pendently. 

Note that special care is needed to prepare a sys- 
tem where a single crystalline cluster coexists with sur- 
rounding fluid with which it is in equilibrium. When we 
would start out with a packing fraction 77 that exceeds r]f 
slightly, such as those values from which data in Fig.[3]are 
included, but choose an initial state that is similar to the 
bulk metastable fluid (without walls) as an initial con- 
dition, the system develops towards a metastable state 
where the density of colloids just show the familiar layer- 
ing at both walls (cf. the analogous data for hard spheres 
at 77 < 77/ presented in [64|), but it would take an unac- 
ceptable long time in the simulation until a wall-attached 
solid cluster would be nucleated. This expectation actu- 
ally is born out by the simulations (see next section). 
This fact is understandable from rough estimates of the 
barriers for heterogeneous nucleation, obtained from the 
standard classical theory [H [TH |T5 for the sizes of solid 
clusters as studied here, which need only the contact an- 
gle (which is estimated independently for our system, 
see below) as an input: according to the classical Turn- 
bull [9l[T0] theory of heterogeneous nucleation the volume 
V* of a critical droplet having a sphere cap shape and 
contact angle (cf. Fig. [4]), the associated free energy 
barrier is 



and hence we can conclude that f{0) = 5/32 (or larger). 
Using then a particle number A^* = 4413 in our solid clus- 
ter (which is a typical example) and taking for the solid 
the packing fraction at coexistence pressure, rjm = 0.64, 
we find that the corresponding volume is F* = A^*/p = 
N*7i(j^ /(Qrjm) ^ 3610 (remember that a = 1 is our 
unit of length). Using then the estimate for f{9) as 
quoted above we find R" ^ 17.67 and AF^^ ^ 20437^^. 
Now the estimate of Zykova-Timan et al. [40 for our 
model {Eq. ^} is 7,^ ^ (0.95 ± Om)kBT/(TK Conse- 
quently, one would predict a free energy barrier as large 
as AFj^g^ ^ 2000/cbT! Even if this estimate would be 
an overestimate by a factor of two or three (which is 
well possible in view of the crudeness of the approxi- 
mations Eqs. (11), (12) for crystal nucleation) still the 



spontaneous formation of such large crystalline clusters 
as studied here would never be visible in a simulation. 

Thus the recipe to study large wall-attached crystalline 
clusters is to prepare the system in an initial state from 
which there is either only a low free energy barrier to be 
crossed for the system on its way toward thermal equilib- 
rium (or, even better, no barrier at all). This is achieved 
by putting a crystalline seed of roughly the right size 
into the box (Fig.[4|. This seed either has the shape of a 
cuboid or hemisphere. In either case an integer number of 
close-packed (111) lattice planes of the fee structure are 
stacked upon each other parallel to the confining wall at 
the bottom of our "container" . The lattice constant of 
this crystalline cluster is chosen such that the crystal has 
a packing fraction rj = r]m- The volume of the crystal is 
cut out from di L x L x D simulation box filled by well 
equilibrated fluid at packing fraction 77/. The lever rule 
then fixes the volume of the crystal nucleus for a given 
choice of r] for r]f < r] < rjm'- 



y* = (47ri?*V3)/(e), AFh^t = AFhU/W (H) 



where 



Vt] = Krystal^m + (V - Crystal )^/ • 



(13) 



^Kom 






R^'^isi 



f{0) = (l-cosi9)^(2+cosi9)/4 . 

(12) 

In Eqs. (11), (12) it is assumed that the dependence 
of the solid-liquid interface tension 7^^ on interface ori- 
entation can be neglected, and it is also assumed that 
75^ does not depend on the radius of curvature R* of 
the "droplet". Then the critical droplet that forms in a 
homogeneous nucleation process has a spherical shape, 
and the free energy barrier against homogeneous nucle- 
ation ^F^Q^ is just (1/3) of the total surface free energy, 
47ri?*^7s/ of the critical droplet [T, ^2]. For heterogeneous 
nucleation, the free energy barrier that needs to be over- 
come, is reduced by the same factor f{0) as the volume 
of the sphere cap is reduced in comparison with the full 
volume of the sphere [H [10] . 

Estimates of the contact angle from Young's equation 
(see below) imply that is at least as large as ^ ~ 70^, 



Of course, Eq. (13) disregards finite size effects: the 



packing fraction of the liquid coexisting with a nanoscop- 
ically small crystal is expected to slightly exceed r]f be- 
cause the density of the fluid surrounding the crystalline 
cluster must be enhanced due to the Laplace pressure as- 
sociated with a small "droplet"; similarly, also the den- 
sity of the nanocrystal may differ somewhat from its 
macroscopic counterpart. However, for the already some- 
what large (A/"* > 1000) nanocrystals studied here, these 
effects did not prevent the successful equilibration of the 
crystalline clusters. As a first step of this equilibration, 
forbidden overlaps of particles in the crystal and in the 
fluid are removed. We found that using a hemisphere 
as an initial state equilibrium is reached more rapidly 
than with a cuboid as initial state; but the properties 
of the equilibrium that is reached do not depend on the 
initial state, as it should be. Of course, in equilibrium 
we could have the solid cluster attached to the wall at 
z = D with the same probability as the situation that 
is actually studied: Due to the initialization (Fig. [4| the 
symmetry between both otherwise identical walls is bro- 



ken "by hand" . Note also that a consideration along the 
lines of Eqs. ( pTj ), ([l3| readily shows that a situation with 
two solid clusters (one cluster at each of the walls at 2: = 
and z = D) is less favorable than the single cluster state. 
In order to study the detailed physical properties of 
the crystalline cluster, it is necessary to identify in mi- 
crostates of the system which particles belong to the 
solid and which particles belong to the liquid. For this 
purpose, we follow the traditional methods [40l [72l [73] 
where the "coherence property" of a particle and its near- 
est neighbors are analyzed using spherical harmonics. 
Specifically, one computes the complex vector qeii) for 
each particle (labeled by index i). The 13 components 
(labeled by m) of this vector depend on the relative orien- 
tation of the "bond" connecting the particle to its neigh- 

bors, and are defined as < Q^mij) = (X/Nhii)) Yl ^em \ 

^ k=l ^ 



q6,m{i) = 



Qemji) 



(E \Q6m{^wy^' 

171 — 6 



(14) 



where Nij{i) is the number of nearest neighbors of par- 
ticle z, k labels the bonds connecting particle i with its 
/c'th neighbor, and Qim are spherical harmonics {1 = 6 
has to be used in the present case). Such nearest neigh- 
bors of the z'th particle are identified by defining a cut- 
off distance, and all the particles whose relative distance 
from the i'th particle is within the cutoff range are iden- 
tified as candidates for being nearest neighbors. The cut- 
off distance is chosen as the first minimum in the radial 
distribution function of the colloid particles. Then we 
compute dQ{i) according to 



deii) 



E 

k—l m= 



+6 



Yl QQ,rn{i) ■ Qemik) 



(15) 



We consider z as a particle belonging to the solid if 
dei'i) > 0.85 and its number of nearest neighbors is 
Nb{^) ^ 12. The latter choice has the consequence that 
the interface between the crystal and the liquid is put 
towards the crystal region of the (extended [40]) liquid- 
crystal interfacial profile, rather than into its center, so 
some roughness from the surface of the crystal is elimi- 
nated, and the number of particles in the crystal slightly 
underestimated. But we expect that this choice will not 
lead to noticeable systematic errors of the contact angle 
of the crystalline cluster. If we choose a smaller value 
than 12 for this cutoff, also small clusters in the liquid, 
which are not of physical significance, are counted as be- 
ing crystalline. 



III. RESULTS AND DISCUSSION 

Figs. [5| [6] show typical snapshot pictures of the crys- 
talline clusters obtained by the method as described 



above (particles identified as fiuid are not shown). One 
can see that we do obtain crystalline clusters of roughly 
sphere cap shape (note, however, that there occur sub- 
stantial fiuctuations in both the size and the shape of 
these clusters, as expected, since we do not apply any 
constraint to the properties of the solid cluster, other 
than that the lever rule, Eq. (13), must be satisfied, since 



the total particle number in the simulation box is a con- 
served quantity). Fig. [t] shows cases where rj was chosen 
too large (for the considered choice of box linear dimen- 
sions) so that no longer a single sphere cap is stable, but 
rather the system forms a (distorted) cylindrical cluster 
(connected in itself by the periodic boundary condition, 
as drawn schematically in Fig. [2| or even a slab-like con- 
figuration forms. 

From Figs. [5]- [T] it is clear, that reliable data on crys- 
talline clusters can only be obtained as long as the lat- 
eral linear dimensions of the crystalline cluster are dis- 
tinctly smaller than the box linear dimensions L^ and 
Ly. We shall disregard conditions where cylinder-like and 
slab-like domains form in the following. The solid clus- 
ters that are not affected by the lateral periodic bound- 
ary conditions are analyzed under the assumption that a 
sphere-cap shape is a reasonable approximation (Fig.[8|. 
Figs. ([9|, ( p^QJ ) show typical results for the time evolution 
of the cluster size A/"*, contact angle 6>, basal radius r and 
cluster height /i, and the resulting probability distribu- 
tions of these quantities. Despite the use of a significant 
computational effort (typically a million Monte Carlo 
steps per particle, for systems containing on the order of 
60000 particles, were used) the amplitude of the fiuctu- 
ations in the crystalline cluster properties are large, and 
the correlation time of these fluctuations typically is of 
the order of 10^ MCS. Thus it is diflicult to ascertain the 
systematic trend that the comparison of Figs. |9] and 10 
suggests, namely that there is a systematic increase of the 
contact angle with the size of the crystalline cluster. Note 
that we have made runs both for a square basal plane of 
the box {Lx = Ly) and for a hexagonal base, but we did 
not find that this choice leads to systematic differences. 
We also emphasize that the total internal energy in the 
system is much less fiuctuating, and its average shows 
a smooth variation with rj (Fig. 11). We have taken all 



precautions that our systems are well equilibrated and 
our runs do constitute a significant statistical effort: but 
the equilibrium between the (sphere-cap shaped) solid 
cluster and its environment, which is only stabilized by 
the constraint of constant density {Eq. ([l3|} and which 
would not be stable in the finite box if one could carry 
out the simulation at constant chemical potential rather 
than at constant density, allows very strong and long- 
lived fiuctuations. There is no constraint on the shape 
of the crystalline cluster other than the driving force to 
minimize the free energy of the total system. The same 
fact holds concerning the size of the cluster: a large fluc- 
tuation increasing the cluster size needs a density fluctu- 
ation in the surrounding fluid which then has a too low 
density. But the driving force to bring the fluid density 



back to its appropriate value is rather weak, and hence 
it may need a long time for the cluster size to return 
to its equilibrium value. Similar considerations apply to 
other cluster properties as well. In any case, these large 
fluctuations of the contact angle, height, basal area, and 
volume of the cluster constitute evidence that the cluster 
surface is rough, rather than faceted: in the latter case 
much less fluctuations would be expected. 




FIG. 5. Examples of typical crystalline clusters, as seen 
from above in a three-dimensional view. The number of solid 
particles in the clusters are AT* =3092 (a), 4743 (b), 2797 (c) 
and 4852 (d). Size of the system used for (a)-(b) Lx — Ly — 
39.60654835212, D = 40.31380814413 and that for (c)-(d) is 
Lx = 45.73370270554, Ly = V3/2L^, D = 40.31380814413. 
Corresponding total volume fractions in the Lx ^ Ly x D 
simulation box are quoted in the figure. All crystals shown 
here have been chosen with their (111) planes parallel to the 
substrate. 
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FIG. 6. Same as Fig. |5] but viewing the projection into the 
xz-plane. 

Given the fact that from the averaging of proper- 
ties of the observed wall-attached crystalline clusters 
(Fig. ^ [To]) one can conclude that the concept of wall- 
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FIG. 7. Same as Fig. |5] but for j] = 0.51539(a) and 

r\ — 0.545(6). In case (a), a (distorted) cylindrical domain 
has formed, while (b) corresponds to a slab- like configu- 
ration. The number of solid particles in the clusters are 
AT* = 8358 (a) and 23823 (b). Size of the system used for 
(a) is Lx ^ Ly ^ 39.60654835212, D = 40.31380814413, 
and that for (b) is Lx 
D = 65.42781457594. 



33.59395945368, Ly = V3/2L^, 



attached sphere-cap shaped "droplets" is at least qual- 
itatively reasonable as a coarse-grained description, it 
makes sense to compare the observed contact angle of the 
"droplets" with the "macroscopic" contact angle. The 
latter results from Young's equation [16-18 as 



Iwf - Iwc = Ifc COS 



(16) 



Here j^f is the excess free energy of the fluid due to 
the confining wall, and j^c the analogous quantity of 
the crystal, for the case that the close-packed planes in 
the fee crystal are parallel to the flat wall surface (con- 
sistent with what is actually observed. Fig. [6]). Note 
that in Eq. (16) it is explicitly assumed that the inter- 



facial tension of the fluid-crystal interface does not de- 
pend on the orientation of this interface, which clearly is 
an approximation, and will not hold true in general. If 
we nevertheless accept this approximation, we can con- 
clude from the work of Zykova-Timan et al, [40 that 
7/c = (0.96 ± 0.05)kBT/a^. For the estimation of the 
wall tensions j^f and j^c recently several fairly accurate 
methods were developed (64] [TOi . 

Fig. 12 plots results for the wall tensions j^f^ Iwc for 
the present model (defined by Eqs. (|6|,- (|9|). Using then 

at the transition 



the resulting estimates A7 = 7. 



w f 7w 



we can use Eq. (16) to predict the contact angle 6>, which 
turns out to be close to 70^. Fig. 12 shows that A7 (and 
hence 0) depend on the strength e of the WCA potential 
only rather weakly. A related finding was already re- 
ported in [64] for the simple hard sphere fluid (for which 
complete wetting seems to occur; i.e. 6> ^ 0^). These re- 
sults imply that the variation of the strength of the inter- 
particle attraction {rf!p) is a suitable recipe to change the 
contact angle of the system, while the variation of the 
strength of the wall repulsion {e) is not. Recall that ex- 
perimentally rfp can be varied by changing the polymer 
concentration in the system, while e could be varied by 



different wall coatings (e.g., using a polymer brush layer 
of variable grafting density) . 
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FIG. 8. Measurement of the contact angle of a crystalline 
cluster assuming a sphere cap shape. When a crystalline clus- 
ter is identified the number of particles Nh in the first crys- 
talline layer adjacent to the basal plane of the substrate is 
estimated. From the area Nhcr'^ of this layer we obtained 
the radius r in the basal plane as r = ^jN^j pAi^ (remember 
cr = 1), here pA is the areal density of the cluster base which 
is equal to pA = ;^(^)^''^- From the total number iV* of 
particles in the crystalline cluster one then obtains the height 
h of the sphere cap via N* — r]mh{3r^ -\- h^)- The contact 
angle then follows as ^ = arccos [(r^ — /i^)/(r^ + /i^)]. 
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FIG. 9. Properties of wall- attached crystal clusters, for a 
packing fraction 77 = 0.508791 of particles in the simulation 
box {L:, = Ly = 39.60654835212 and D = 40.31380814413). 
Left panel shows the time variation of the cluster size (A/"*), 
contact angle (0), basal radius (r) and height (h). The right 
panel shows the resulting distribution functions of these quan- 
tities. 

Fig. [13] then summarizes our results for wall-attached 
crystal clusters, plotting contact angle ^, basal radius r, 
particle number of the cluster A/"* and height h versus 
the chosen total packing fraction 77 in the system. Two 
different box sizes are used: as expected, A^* (and also 
h and r) must increase when the box size increases: as 
expected from Eq. (13), taking the box size to infinity 
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FIG. 10. Same as Fig. ^ but for r] = 0.507497 (L,, 
45.73370270554, Ly = V3/2La: and D = 40.31380814413). 
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FIG. 11. (a) Time evolution of the total internal energy per 
particle for several choices of packing fraction 77, as indicated 
and (b) the average total internal energy per particle plotted 
vs. packing fraction. Size of the system used for all the choices 



of packing fraction is Lx 
40.31380814413. 



45.26462668814 and D 



at constant 77 in the two-phase coexistence region we ap- 



proach macroscopic phase coexistence. Approaching this 
limit, the contact angle 6 should not change. Gratify- 
ingly, we do find that the two data sets yield contact an- 
gles ^ = 70 ± 2 degrees irrespective of the choice of box 
size and packing fraction. Of course, for very small crys- 
talline clusters (containing a few hundred particles only) 
a systematic effect on the contact angle is expected due 
to line tension effects, for the cluster sizes shown (where 
A^* is several thousands) such effects are too small to be 
distinguished, given our statistical errors. Of course, in 
view of all the possible uncertainties whether or not jfc 
depends on the interface orientation, and consequently 
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FIG. 12. Wall tension of the AO model plotted vs. bulk 
colloid packing fraction 775, as obtained from the "ensemble 
mixing method" , as described in [BU [70] . Two choices of 
WCA potential are included, and the resulting differences 
A7 = jwf — Jwc at the phase transition are quoted. Note 
that the horizontal straight lines indicate the two-phase co- 
existence region. 



whether a sphere cap shape is accurate, etc, the agree- 
ment of this finding for with the corresponding pre- 
diction based on the Young equation may be somewhat 
accidental. More work on this problem is desirable. 

In the studies of phase coexistence at the vapor-liquid 
transition [Ml HSl US] |46], it was shown that it is also 
useful to analyze droplet properties as a function of the 
chemical potential of the vapor surrounding the liquid 
droplet, since then a part of the data collapses on "mas- 
ter curves" that are independent of the box linear di- 
mensions (the parts that do not collapse are still af- 
fected by undesirable finite size effects due to the droplet 
evaporation-condensation transition for small droplets or 
by the transition from spherical to cylindrical shape for 
large droplets). In the present case of a liquid-solid tran- 
sition, due to the high density of the liquid, the chemical 
potential of the liquid is not straightforwardly estimated, 
but we can easily determine the packing fraction rjb of 
the bulk metastable liquid that coexists with the crys- 
tal cluster (this is done by sampling from a slab near 
z = D/2). 

Fig.fg 



Using the data of Fig. 13, this is done in 



We expect N* (as well as r and h) to decrease 
monotonically with increasing 775 (the smaller the clus- 
ter, the denser the surrounding fluid must become, due 
to its pressure increase by the Laplace pressure which 
increases proportional to the inverse radius of curvature 
of the cluster). We see that most of the data do fol- 
low this expectation, but part of the data for the smaller 
system (for A/"* > 7500) break off from the common, size- 
independent "master curve" and also the corresponding 
contact angle data get systematically smaller. Similar 
trends were also seen when too large droplets were in- 
cluded into the analysis of the vapor to liquid transition 
[m [151 [45l [46] : iii this case it could be clearly proven that 
this trend is due to the problem that occasionally the 
droplet undergoes a transition in its shape from spheri- 
cal to cylindrical [45l [46] (or sphere-cap to cylinder-cap 



in the presence of walls [Ml [15] , respectively) . This prob- 
lem is easily missed without careful analysis of time evo- 
lutions of droplet properties, such as shown in Fig. [9J 10 
(which refer to smaller A/"*, where this effect did not yet 
occur). Similarly, also the data for the larger box size 
show a dramatic rise for 7^5 < 0.5035, where TV* > 10000, 
and we feel that these data should be discarded for the 



same reason. Of course, this analysis of Fig. 14 is a first 
step only: it clearly would be desirable to extend this 
study to both smaller and larger box sizes, but in view 
of the huge demand in computer resources needed, this 
is left to future work. When we tentatively discard these 
data where N* is presumably too large, we obtain some 
evidence that decreases with increasing 7^5 (and hence 
decreasing r) slightly. 
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FIG. 13. Properties of wall-attached crystalline clusters, 
shown as functions of the packing fraction ry: contact an- 
gle (a), basal radius r (b), particle number N* in the 
crystalline cluster (c) and height h (d). Two choices of to- 
tal box size are included: the smaller system has Lx — 
Ly = 45.26462668814, /:> = 40.31380814413 (red color data 
sets with O symbol), the larger system has Lx = Ly = 



50.92270502416, D 
with n symbol). 
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FIG. 14. Same as Fig. |13| but shown as functions of the 
packing fraction r]b of the bulk metastable liquid coexisting 
with the crystalline clusters. 

From the knowledge of the contact angle and of the 
radius i?* of a droplet one can make a prediction of the 
corresponding free energy barrier AF^^^ that needs to be 
crossed, if such a droplet forms spontaneously by thermal 
fluctuations. 
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FIG. 15. Schematic description of the loop in the pressure 
versus packing fraction isotherm at a hquid-to-sohd transi- 
tion in a finite simulation box. In the thermodynamic limit, 
p rises with 77 up to the value pco at 77 = 77/, and stays con- 
stant at p = Pco throughout the two-phase coexistence region, 
until it continues to rise at 77 = ?7m- In a finite system, the 
region of homogeneous fluid is enhanced up to ry = 77^ , where 
the droplet evaporation/condensation transition occurs, and 
then it decreases smoothly, until another transition occurs 
from spherical to cylindrical shape of the droplet. (In sys- 
tems with walls, the actual droplet shapes are sphere-cap like 
or cylinder-cap like, respectively). Only for the region of 77 
where the two-phase coexistence in the finite box corresponds 
to a crystal slab separated by flat walls from the fluid is the 
pressure (at least approximately) equal to pco- For larger 77 
the roles of fluid and crystal are interchanged. Note that in 
reality the transitions in Fig. [15] are not sharp but slightly 
rounded. When the volume of the system tends to infinity, 
T/t -^ r/m, and the pressure enhancement p{r]t) — Pco -^ 0, 
so while the extrema of the loop in Fig. [15] gets the sharper 
the larger the system becomes, the whole loop ultimately has 
disappeared after the thermodynamic limit has been taken. 



However, experience with nucleation in vapor-liquid 
transitions or fluid-fluid unmixing [Ml |l5l |45-47 sug- 
gests that barriers predicted from Eqs. (11), (12) often 



are significantly larger than the actual barriers. In order 
to test whether this problem also occurs in the present 
case, we recall that classical nucleation theory predicts 
also a relation for the pressure p' of the liquid coexisting 
with the droplet, namely [U H] (cf. Fig. flS] for notation; 
Fig. 15 is a schematic counterpart of Fig. Isj) 



p'-p,o = {2^fc)/R* 



(17) 



Using our estimate for i?* and 7/c from [40] one 
predicts that for i?* = 14.4 the pressure difference is 
p' — Pco ^ 0.13 while the actual pressure difference seen 
for this case in Fig. [sjis p' — Pco ^ 0.4. Potential causes 
for this failure include the assumed spherical shape of 
the crystalline cluster in eq. [17] and the exact definition 
of R* in our simulation. Clearly, more work is required 
to study the surface free energy of the curved crystalline 
clusters. 



IV. CONCLUSIONS 

In this paper we have presented a study of liquid-solid 
phase coexistence in the constant volume (NVT) ensem- 
ble for the Asakura-Oosawa model of a colloid-polymer 
mixtures, focusing on the case of the size ratio q = 0.15 
and polymer reservoir packing fraction rj^ = 0.1, for 
which the bulk phases and the interfacial stiffness 7100 
for a fluid-solid interface with an (100) surface of the 
crystal were studied in previous work [40]. Using rather 
large systems (containing of the order of almost 10^ col- 
loidal particles), we have shown that an analysis based 
on lever-rule type arguments allows us to gain insight on 
many aspects of phase coexistence, both with respect to 
bulk properties characterizing it, and with respect to the 
contact angle of sphere-cap shaped crystalline clusters. 
When we choose the packing fraction rj roughly half way 
in between the bulk coexisting liquid {r]f) and solid (rim) 
phases, we obtain a phase coexistence, where adjacent 
to the left wall there is a crystal film (with (111) planes 
stacked on top of each other up to a distance z = Di), 
followed by a liquid up to the distance D — D^^ and then 
another crystal film up to the right wall (at distance D 
from the left wall) follows. We have checked (Appendix 
A) that this liquid slab that forms in between these crys- 
tal layers is in full thermal equilibrium and its packing 
fraction 77^ and pressure Pco can be measured very pre- 
cisely in this way. One can show that these values are 
independent of r] (and hence direct evidence for the flat 
horizontal portion of the pressure vs. -packing fraction 
isotherm. Fig. [sj is provided), and from the observation 
of the profile r]Xz) , Fig. [im one can also measure Di + D^ 
and thus establish that the lever rule holds, as it should 
be. From the profiles one can also extract the distance d 
between the crystalline (111) planes and thus check the 
self-consistency of the estimation of rjm- 

The main interest of this paper, was the study of wall- 
attached crystal clusters (Figs. [5| [6| [9] [lT| [13] [14] ), which 
were created by special choice of initial states (Fig. [4| 
for packing fractions r] that exceed rjf only slightly. If 
this excess is too small, such crystal clusters were found 
to be unstable and dissolve again, and one is left with a 
somewhat compressed uniform fluid (apart from the lay- 
ering at the walls, similar to what is seen in Figs. [TtJ [18] 
at the right wall). If this excess is too large, crystalline 
clusters of cylindrical shape (or even planar crystalline 
films) form (Fig. [t]). Using the liquid- wall and crystal- 
wall surface tensions that were estimated by a different 
method [64 and are shown in Fig. 12, one can estimate 
the contact angle to be close to 70^ (this estimate relies 
on the assumption that the fluid-crystal interface tension 
7/c is approximately independent of crystal surface orien- 
tation, and hence the estimate of [40 can be used here). 
Gratifyingly, the numerical data observed for the con- 
tact angle (Figs.[9j[l0j[l3J[l4| are compatible with this 
estimate. Of course, the crystal-wall interfacial tension 
7^(;c is expected to depend significantly on the orientation 
of the crystal axes relative to the wall. Consistent with 
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this expectation, a contact angle close to 90^ is observed 
(Fig. [l7|) if (100) planes of the crystal are chosen to be 
parallel to the wall surface. The estimate of ^y^c foi" the 
(100) orientation from the "ensemble mixing" method is 
live — 1.82 ± 0.051 which leads to a contact angle of 

e ^ 90^ 

We have also pointed out that in the finite system the 
pressure versus packing fraction isotherm (Figs.lSJIlS]) ex- 
hibits a loop, which has nothing whatsoever to do with 
van der Waals-like loops; however: it is entirely caused by 
interfacial contributions to the free energy of the system, 
and since the latter are down by a surface to volume- 
ratio in comparison with the bulk, the loop gradually 
develops towards a fiat variation p = Pco from r] = r]f to 
V — Vm^ when the thermodynamic limit is taken. Since 
the actual pressure enhancement (in the region where p 
decreases with increasing 77, Figs, pj fT3| should contain 
information on the interfacial free energies of the (spher- 
ical or cylindrical crystalline clusters), we have tried to 
study this enhancement of the pressure, too, but using 
Eq. (17) did not yield results that were quantitatively 
consistent with our other results. This problem clearly 
requires further study. 

We like to emphasize, that the present work is a first 
step only; it is necessary to develop thermodynamic 
integration-based methods, from which the excess free 
energy of the system (due to the crystalline cluster) can 
be reliably extracted, in order to be able to make quanti- 
tative predictions for nucleation barriers in such systems. 
Also it will be interesting to study the same model for 
other values of q and r]'^. In any case, it would also be 
very interesting if experiments on colloid-polymer mix- 
tures, which are rather well described by this model, were 
performed. 
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FIG. 16. Plot of the packing fraction profile r]{z), ordinate 
labels on the right side of the picture, and of the local es- 
timates d{z) for the distance between neighboring planes in 
the fee ABCAB- • • stacking, ordinate labels on the left side of 
the figure, for the choice of linear dimensions Lx = 33.59396, 
Ly = 29.09322, D = 64.83302, ry = 0.55(a) and D = 63.67528, 
rj = 0.56(b). The thick vertical lines in the center of these fig- 
ures mark the part of the liquid slab that has been used for 
the estimation of 77/ and pco ■ Note that this figure shows only 
two examples out of many more that were recorded. 



Appendix A: Phase coexistence between 

wall-attached crystal films with a liquid slab in 

between 



When we prepare a system at a packing fraction about 
halfway in between the packing fractions 77/, rj^n of the 
bulk coexisting phases, we expect that the final equilib- 
rium state will be a layered structure, with a crystal film 
of thickness Di attached to the left wall, a crystal film of 
thickness Dr attached to the right wall, and a fiuid slab 
of thickness D — D^ — Di in between. From symmetry, 
one expects D^ = Di^ of course, and the total thickness 
of the crystal films Dr + Di is fixed by the lever rule. 



Such a state is easily obtained if we first prepare the 
total Lx X Ly X D system as a crystal of packing frac- 
tion rjrn (choosing linear dimensions Lx, Ly, D such that 
any misfit with the fee crystal structure at r]m is strictly 
avoided), and removing then particles from the region 
near z = D/2 from the system, until the desired average 
packing fraction is obtained. Then the system is equili- 
brated, and one observes that quickly a liquid slab forms 
in between the crystal layers (Fig. 16). The thicknesses 
of the two crystal layers Dr, Di are only roughly equal 
to each other, but as long as each crystal contains many 
(111) layers stacked parallel to the wall, no noticeable 
systematic error is caused by this slight asymmetry. One 
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sees that irrespective of the precise value of 77 that is cho- 
sen, one does obtain an extended region near z = D/2 
where the volume fraction profile r]{z) is flat, and this 
horizontal region yields an accurate estimate for rjf. Note 
that for a wide range of choices of 7^(0.53 < r] < 0.58) and 
also for several choices of the linear dimensions we always 
obtain the same value of r]f ~ 0.498 (slightly larger than 
the estimate r]f = 0.494 of Zykova-Timan et al. [40 ). 
However, the coexistence pressure Pco (extracted from 



the region when r]{z) in Fig. 16 is flat, using the method 
described in [70]) agrees with the previous estimation [40] 
Pco = 8.00 ±0.01 within the statistical errors (the present 
estimate is Pco = 7.98±0.01). Starting from the initial es- 
timate for the packing fraction of the crystal rjm = 0.64, 
one obtains for the distance between the crystal (111) 
planes d = 0.857. Within the accuracy with which d 
can be estimated from the profiles r]{z) in the crystalline 
films, this is the value that the simulation yields, Fig.[TTj 
Of course, one cannot reliably measure d right at the 
walls, and one should also avoid using the profiles in the 
region of the crystal-liquid interface. Estimating Di -\-Dr 
from the location of the interface positions in the profiles 
r]{z)^ we have confirmed Eq. (11) quantitatively. 

This type of crystal-liquid coexistence simulation 
hence not only provides direct evidence for the strictly 
horizontal part of the pressure versus packing-fraction 
isotherm in Fig. flSJ but yields direct estimates for 7^/, 
r]m and Pco with very good precision. Of course, we have 
simplified matters by using the (previously known) value 
of rjrn to choose linear dimensions L^ , Ly commensurate 
with the fee lattice structure that the system wants to de- 
velop. If the initial choice of 77^ would be somewhat off, 
one would find that the crystal structure exhibits some 
elastic distortion: the distance between planes would 
come out either somewhat larger or smaller than pre- 
dicted from the (wrongly chosen) initial value for rjm'- 
then an iterative improvement of this choice would be 
necessary. Thus, we propose such studies of phase coex- 
istence as an additional method to precisely characterize 
liquid-solid transitions. 



Appendix B: Simulation of metastable crystals with 
(100) planes oriented parallel to the walls 

In the main text, we have described the technique to 
prepare crystalline clusters where by construction of the 
initial states the familiar ABC ABC stacking of the close- 
packed (111) planes parallel to the planar walls ware 
created. It is possible, however, to choose initial states 
where instead the (100) planes are stacked parallel to 
the planar walls. While in the (111) planes the particles 
form a triangular lattice, where (for 77^ = 0.64) the lat- 
tice spacing is 1.0497, and the distance between planes 
is d = 0.857, for the (100) planes the particles form a 
square lattice, with the same lattice spacing, but the dis- 
tance between the planes is slightly smaller, than for the 
(111) stacking, namely d = 0.7423. Figs. 17, 18 compare 
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FIG. 17. (a) Packing fraction profile and (b)-(c) two different 
perspectives of a crystalline cluster with (100) plane stacking 
for J] = 0.513347 



typical cases of clusters with (100) stacking and (111) 
stacking. Note that at the right wall (z close to D^ there 
is the typical layering of the liquid phase near a flat re- 
pulsive wall, which is very similar in both cases. In the 
profile near the left wall, the first peak of r](z) adjacent 
to the wall is again in part due to the layering of the 
fluid and in part due to the crystalline cluster, and again 
similar in both cases. However, while for (111) stacking 
the further density oscillations (which are mostly due to 
the crystalline cluster) decrease monotonically with the 
distance z from the wall, this is not the case for (100) 
stacking: the 3^^ and 4*^ peak of the oscillations are 
less high than the 5^^ to 8^^ peaks. This non-monotonic 
behavior of the peak heights has an obvious interpreta- 
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tion in terms of the cluster shapes: for (100) stacking 
the contact angle exceeds 90^ slightly, and so the cross- 
sectional area of the crystalline cluster along the 5*^ to 
8*^ plane is slightly larger than along the 3^^ and 4*^ 
plane. The bottom snapshot gives direct visual evidence 
for this interpretation. Due to the larger contact angle 
this crystalline cluster is only metastable. 




Ti = 0.51 2462 



(c) ^^^ 



FIG. 18. Same as Fig. 17 but for (111) stacking and 77 
0.512462 
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